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Abstract. We present numerical renormalization group (NRG) calculations for a 
single-impurity Anderson model with a linear coupling to a local phonon mode. We 
calculate dynamical response functions, spectral densities, dynamic charge and spin 
susceptibilities. Being non-perturbative, the NRG is applicable for all parameter 
regimes. Our calculations cover both weak and strong electron-phonon coupling for 
zero and finite electron-electron interaction. We interpret the high- and low-energy 
features and compare our results to atomic limit calculations and perturbation theory. 
In certain restricted parameter regimes for strong electron-phonon coupling, a soft 
phonon mode develops inducing a very narrow resonance at the Fermi level. 



1. Introduction 

The effects of lattice vibrations on the electronic properties of metals is a well studied 
branch of condensed matter physics. These effects include specific heat enhancements, 
strong thermal contributions to the electrical resistivity, and in some materials, 
superconductivity at low temperatures. There remain, however, some outstanding 
problems in the field. For example, the extensive studies of small polaron formation, 
based on the model from the pioneering work of Holstein |l|], have given a very precise 
understanding of systems with one or two electrons present, but there are no comparable 
results for systems with a finite electron density. There are also many strongly correlated 
electron systems, such as heavy fermions or valence fluctuations systems, where there is 
known to be strong coupling between the electrons and the lattice vibrations, in addition 
to the strong inter-electron interactions. The interplay of phonon and electron effects 
in these systems has not received much attention. The reason for the lack of progress 
with these problems is due to the fact there are no reliable calculational techniques 
in the strong coupling regime, where standard perturbational methods break down. 
Some non-perturbational techniques have been developed and successfully applied to 
strongly correlated impurity models but in general these cannot be extended to models 
of lattice systems. However, with the introduction of djTiamical mean field theory 
(DMFT) 01, it has proved possible to extend some of these non-perturbative techniques 
to a wide class of lattice models. This approach exploits the fact that certain infinite 
dimensional lattice models can be mapped onto effective impurity models, together with 
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a self-consistency condition [0]. The calculations for the effective impurity models can 
be carried out using an appropriate non-perturbative method and iterated until self- 
consistency is achieved. Though the method is exact in the large dimensionality limit 
for an appropriately scaled model , the method can be applied to a much wider class 
of models and there is evidence that it provides a very good first order approximation 
for three dimensional systems. (It clearly breaks down for one- dimensional systems but 
one-dimensional systems are a class apart.) 

Not all techniques for solving strongly correlated impurity models can be used for 
the effective impurity models generated through the application of DMFT. The methods 
must be able to handle the dynamical properties of the impurity models, which rules out 
techniques such as the Bethe ansatz. The method which has been used most extensively 
for this type of calculation is the quantum Monte Carlo method (QMC). This approach, 
however, does have the drawback of not being capable of being applied to the very 
low temperature regime. A method which can be used at T = is the numerical 
renormalization group approach (NRG). As developed originally by Wilson for the 
Kondo model [^, |^, ^ it was used only for the calculation of thermodynamic quantities 
but the approach was extended later to the calculation of the dynamic properties of 
impurity models 0, |^, which is required to be able to use the method in conjunction with 
DMFT for lattice models. The DMFT approach together with the NRG has been applied 
to the half- filled Hubbard model to study the metal-insulator transition at T = P, p!0|]. 



and these calculations have recently been extended to finite temperatures |TT|. Our goal 
is to extend this approach to include the effects of coupling to local lattice modes, and 
address some of the outstanding questions concerning the interplay of strong correlation 
and strong electron-phonon interactions. The results should provide some insights into 
the physics of systems such as heavy fermions, valence fluctuators, manganites, and the 
cuprate high temperature superconductors. 

In this paper we perform NRG calculations solely for an impurity model with a 
coupling to a local phonon mode. This is a necessary first step to the eventual aim 
of incorporating the self-consistency constraint of the DMFT and the application to 
models for lattice systems. We introduce the impurity model in the next section. 

2. The Anderson-Holstein Impurity Model 



The model we study is a single impurity Anderson model |12| with a linear coupling to 
a local phonon mode, as in the the Holstein model We will refer to this model as 
the Anderson-Holstein (A-H) model. The corresponding Hamiltonian for the model is 

H=Y: efflf^ + UflfJlU + A(6^ + 6)(E flf. - 1) 

k,(T kcr 

which describes an impurity d level e^, hybridized with conduction electrons of the host 
metal via a matrix element \4, with an interaction term U between the electrons in the 
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localized f (or d) state, where n^o- = f^f^ is the occupation number of this state which 
is linearly coupled to a local oscillator of frequency ujq. 

A very similar model which does not have the interaction term U, but does not 
quite correspond to the A-H model with U = as the electrons are spinless, has been 
studied previously |T^. Both the model with spin and the spin-less model can be 
solved exactly in the case of one electron only in the system |P!3|. As the interaction 
U does not come into play in this situation the results of both models are essentially 
equivalent. However, there is no exact solution for the many-electron case even for the 
spin-less model and U = 0. 

We will solve the model (|l|) using the numerical renormalization group method. We 
will briefly outline the general NRG scheme which is applicable to any impurity model 
with local interactions and coupling to local bosonic or fermionic fields [|, P]. The 
conduction band of the model is logarithmically discretized and mapped onto a semi- 
infinite linear chain coupled to the impurity at one end such that the hopping matrix 
elements along the chain fall off exponentially (~ A~"/^). The many-body states and 
energy levels are calculated iteratively starting from the impurity adding an extra site of 
the chain in each step. After about 5-6 steps, the number of states has to be truncated to 
the 300-800 lowest ones depending on the application in each step. This method has been 
generalized to the calculation of dynamic as well as static properties 0) §• In applying 
this approach to the model with phonons (Eq. (|1|)), a truncation of higher-energy steps 
is already necessary at the initial step. The probability distribution function of finding 
n excited phonons in the system has a clear maximum around n = n = v/uq, a.s can be 
read off equation and falls off rapidly for larger values of n. So it is sufficient to 
keep only states with n < Ucutos where we have chosen an initial ncutofr = 4n. We have 
used the discretization parameter 1.7 < A < 2., and kept at least 600 states in each 
iteration (see Ref. |^ Ol)- 

To calculate the spectral function on the real axis, we need the Fourier transform 
of the double-time one-electron Green's function for the local d-electron. We use this in 
the standard form, 

GA^) = ^\ ^ , y (2) 

where we have taken the hybridization function A(co') = 7!'J2k\Vk\'^S{u! — e^) to be 
independent of u, corresponding to wide conduction band with a flat density of states. 
The self-energy So-(ti;) will be calculated, as in reference from the ratio of two 
Green's functions. Due to the additional interaction term with the phonon there are 
two contributions so that So-(c<j) = S^(co') + S^(co'), with 

where Fo-(a;) and Mcr{uj) are the higher order Green's functions, 

F^u) = {{fJtf,, fl))^, M^uj) = {{f^b^ + b), fl))^. (4) 
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Details of the derivation of this result from the equations of motion are given in the 
Appendix A. In Appendix B, more details regarding the calculation of M„{uj) are 
presented. 

We will also calculate the spectral density of the Green's function of the phonon 
mode, D{u) = {{{b'^ + b), {b^ + b)))^. This can be expressed in terms of the local charge 
susceptibility, 

D{uj) = D\uj) + \^D\uj){{d,d))^D\uj), (5) 

where 

is the non-interacting phonon propagator, and ((0,0))^^ with O = J2cr fafa — 1 is the 
local dynamic charge susceptibility. 

It will be useful in interpreting the numerical results to compare them with 
results from perturbation theory in the relevant regimes. We will set up the general 
perturbational approach in this section and give specific results later. Starting from 
expression for the partition function Z in the form, 

Z/Zo=(re-/o''^-W'^-)o, (7) 

where T is this usual time-ordering operator and H[^t{T) is the interaction part of the 
Hamiltonian (including the hybridization term) in the interaction representation with 
T as the imaginary time variable and f3 = 1/T, and ( )o denotes a thermal expectation 
value with respect to the non-interacting part of the Hamiltonian. We can remove 
explicitly the phonon operators and rewrite (|^) as 

Z/Zo = (Te-/.f(^-(^)+^-W)^-)o, (8) 

where H^^ is the purely fermionic interaction part of the Hamiltonian ([^) and i^ret is 
an additional retarded local fermionic interaction due to the elimination of the phonon 
terms, which is given by 

Hretir) = \' r d{T)d{r')D\T - r')dT', (9) 
Jo 

where D^[t) is the non- interacting phonon Green's function in the imaginary time form. 
Details of this derivation are given in Appendix C. 

It will be convenient to convert this expression to a path integral over fermionic 
Grassmann variables and integrate out the conduction electrons. We can then write (|^) 
in the form, 

Z = J U'^iU)ViU)e-fo'^^^(^)^^, (10) 



where 



Aff(r) = / rfr'5:/.(r)[G(°)(r - r')]-'Mr') + Un^{r)n^{T) 
+ ^ ['dr'D\r - r'){Y.n.{r) - 1)(E^.'(^') " 1 
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where faij), faij) are Grassmann variables, n^ij) = /o-(t)/o-(t), and G'^^\t) is the 
imaginary-time of the non-interacting one electron Green's function for the localized 
f-electrons. By introducing a Grassmann source field for the f-electrons one can convert 
this partition function into a generating function for the interacting local Green's 
function (jo-(r) in the standard way []T6|. This gives a convenient way of generating 



the perturbation theory and Feynman diagrams for this Green's function. 

In the limit uj q ^ oo such that v = X'^/ujq remains finite the effective interaction 
becomes an instantaneous one equal to —v{J2a''^cr — 1)^- This interaction term can be 
absorbed into the standard Anderson model with the terms linear in n„ being absorbed 
by modifying e/ — > ^f+v, and the remaining part by changing U to U—2v. In this way of 
taking the limit ujq — * oo, the ratio v/ujq —>■ 0. But there is also a high frequency regime 
of interest in which both and v are large such that v/uq ~ 1. For this parameter 
range it is convenient to apply the canonical transformation U~^HU with If given by 

which is a displaced oscillator transformation such that 

b = U-'bU = b-—{Y.n,-l), and = If-^ f^U = e-^^"'-'^ f,. (13) 

^0 a 



H' = U-'HU = Y: tfflf, + Uflfjlf, 

<J 

k,(T ^ ko- 



(14) 



where e/ = e/+f and U = U — 2v. The exponential term in the phonon operators, which 
multiplies the hybridization term, arises because a cloud of phonons may be created or 
absorbed when the f-electron occupation changes, due to the coupling to the lattice. If 
we take the expectation value of this Hamiltonian in the zero phonon state |0) then, as 
^Q|g±A/a;o(6t-b)|o) = g-f/2<^o^ i-j^g effective resonance width becomes A = e'^'/'^oA. For 
f/ciJo ~ 2 or 3, this would amount to a significant narrowing of the resonance for the 
non-interacting model. This exponential renormalization factor is a common feature of 
strongly coupled electron-phonon problems but, as was shown in earlier work of Hewson 
and Newns [ll3| , for the spinless model, this resonance narrowing does not in general 
occur in the strong coupling regime for this model, but only in a certain parameter 
regime. 

In the zero hybridization limit for the transformed Hamiltonian in equation (0) 
the electrons and phonons become decoupled. However in this limit the calculation of 
the Green's function in terms of the original f-electron operators is not trivial, as it 
involves the unitary transformation U but it can be evaluated exactly by using the well 
known techniques for calculating the expectation values of exponentials with a linear 
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combinations of bose creation and annihilation operators [0 . The result for T = is 

^0 ^- ^ uj-lf-nujQ u -ef - U -nujo ^^^^ 

u — if + nujQ u — If — U + nujQ _ 

where the expectation values, such as ((1 — ?2/,o-)(l — are calculated using (0) 

in the 'atomic' limit = for T = 0. 

In the weak coupling regime we can use perturbation theory based on the original 
form of the model given in equation (|lD. Using the finite temperature perturbation 
expansion for the self-energy E^(a;), carried out to order A^, and analytically continued 
to real frequency cu, we obtain 

9 X2 

J [u; — e + ti^o uj — e — ujq J 

where nb(e) and /(e) are the Bose- Einstein and Fermi-Dirac distribution functions. The 
occupation value n/^o- of the f-state with spin cr, and the corresponding spectral density 
Pf^fj{uj), can in principle include self-energy contributions from the local interaction term 
U . These quantities to zero order in U are given by 

Substituting these expressions into (ITB) for the self-energy at T = 0, we find 



TT y \\uj + ujq\ ) uj — ef + ujq — \\uj — ujq\ J uj — Ef — ujq — ] 

and 



(18) 



S^(2)(c.) = + I ~Pf.m,-Ae)Pf,-Ae') ^^''"^^f^ J ede'de", (20) 

J ijj — e + e — e 



^ (cu + cuo - e/)2 + A2 ^ -^(^_^Q + ey)2 + 

The corresponding expression for S^(ci;) to order f/^ is 

^(e,e^e-- 
' - e + e' - e" " 

where 

/^(6, e', e") = /(eO - /(6')/(e") + /(e)/(e") - /(e)/(e'). (21) 

The effects of phonon scattering could be included in the f-spectral density pf^^i^), and 
Uf^a, in evaluating (^). If terms of order were to be included they would generate 
corrections of order A^f/^. These contributions would be distinct from those of the same 
order generated by including terms of order U"^ in evaluating S'^(co') from equation ([T6|). 
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To zero order in A we use (|T^) in evaluating (pOl). These results can be used to compare 
with the NRG ones in the weak coupling regime. 

There is a large parameter space for this model so we will present the NRG results 
in four sections, corresponding to the symmetric and asymmetric models with U = 0, 
and then the same models with U 0. 



3. Symmetric ^7 = Model 

In the upper panel of figure we present the NRG results for the spectral density of the 
f-Green's function for various values of A. The bare phonon frequency has been taken to 
be = 0.05 in this and all subsequent plots. The significant features are the narrowing 
and virtual collapse of the central resonance with increasing coupling strength A and the 
development and outward shift of high energy peaks. This shift is proportional to A^ and 
can be understood by looking at the corresponding results in the atomic limit (V = 0). 
The spectral weights and positions of the atomic-limit excitation peaks of equation (p!5|) 
for selected values of A are plotted in the middle panel of figure |I[ The high-energy 
features in the upper panel correspond to a broadening those seen in the atomic limit. 
In this case (Tf = v and U = —2v, so the effective interaction is an attractive one. 
The peaks below the Fermi level correspond to removing an electron from a double- 
occupied site, the peaks above to adding an electron to a previously unoccupied site. 
The probability that the impurity is singly-occupied for finite A vanishes exponentially 
with A^ (see equation (|TBp). The narrowing of the central peak, the zero-phonon peak 
occurs very rapidly with increase of A^. On the basis of equation (p!^, one might 
expect the narrowing to follow the exponential form ~ e~'"^^° . In figure ^ we plot the 
relative width of the central peak as function of A. Only for small values of A is it 
approximately proportional to exp{—v/uJo). It decreases much more rapidly with higher 
coupling strengths. This could be due to the effect of the softening of the phonon mode. 

The imaginary part of the phonon Green's function, displayed in figure ^, shows 
clear evidence of the emergence of a soft phonon mode already for intermediate coupling 
strength. The inset shows the same as the large figure, but over a larger frequency range, 
displaying also the 'bare phonon' peaks at u = ±0.05. With increasing A, another peak 
appears at lower frequency. This soft phonon mode can be traced back to features in 
the charge susceptibility. Equation (j^) shows the relation between the phonon Green's 
function and the dynamic charge susceptibility. The soft phonon mode reflects the lower 
energy peak in the charge susceptibility. Figure ^ shows the position of the lower energy 
peak in phonon propagator. It falls off approximately with tux = 1/30 exp(—-^^^^^^^). 
With this result, we also reach a good fit to the data of figure ^ if we assume the 
zero-phonon peak width to behave like exp(— f + instead of exp(— f /c^o)- 

A comparison of our results with those of lowest-order perturbation theory, 
equation (IB), which is displayed in the lowest panel of figure reveals the limitations 



of the perturbational result in this parameter regime. We would not expect to get 
the exponential narrowing of the zero-phonon peak from a weak-coupling result, but 
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Figure 1. Spectral density A{uj) = -i9G^(w) for J7 0, ej = 0, F = 0.1 ^ A w 
0.016 and ojq = 0.05. The respective values of A are given in the graph. The upper 
panel shows NRG results, the two lower panels the atomic limit and the perturbation 
theory (PT) rcsidts. For the atomic limit calculation, the height of the peaks is a 
measure of the respective spectral weight of the excitation. 



the high-energy features are also not well reproduced. Although there are high-energy 
peaks, their positions are in poor agreement with the NRG results and their shift to 
higher frequencies is proportional to A instead of A^. 
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Figure 2. Width of central resonance ('zero-phonon peak') as function of A for the 
parameters of figure |l] (NRG calculation) . The solid and dashed lines are explained in 
the text. 
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Figure 3. Imaginary part of the phonon propagator for the same parameters as used 
in figure |l|. The inset shows the same quantity on a larger energy scale. 



4. Asymmetric U = Model 



As soon as we introduce an asymmetry into the model, we find a qualitatively different 
picture. In figure |^ we plot the NRG, atomic limit and perturbation theory results for 
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the spectral density of the local Green's function for a range of values of A similar to 
figure |1|, but with an /-level located at e/ = —0.05. In this spectrum there is a single 
peak corresponding to an excitation from a doubly-occupied impurity site, positioned 
at e/ + f7 — nuo = ej — v — nujQ, where n = v/ujq is approximately the mean number 
of phonons excited. The shift in the peak to lower values reflects that seen in the 
atomic limit, equation (jl^l). Additionally, the peak broadens with increasing A due to 
the higher probability of exciting multiple phonons on removing the local electron. As 
the peak moves away from the Fermi level, the charge fluctuations are suppressed. This 
is manifest in the absence of the soft phonon mode, as can be seen in figure ^ showing 
the phonon propagator in this parameter range (cf. equation (||)). 

In contrast to the symmetric case, the perturbational result, shown in the lower 
panel of figure |^, does give quite good agreement of the peak positions, the increase of 
broadening, however, is absent. 

5. Symmetric f/ 7^ Model 

In figure ^, we plot the spectral density for U = 0.1 and e/ = —0.05. In contrast to 
the U = case, the central peak initially broadens with increase of A and then narrows 
rapidly. When A = 0, U/nA = 2.03, we are initially in the Kondo regime and so we 
can interpret the central peak here as Kondo resonance with a width of the order of 
2Tk ~ {UA)2 exp(^). Hence if U is replaced hj U = U — 2v a.s in equation (0) 
then Tk would be expected to increase, provided U remains sufficiently large so that we 
stay in the Kondo regime. The situation will change as tJ decreases below ttA. When 
U becomes negative we expect the results to be similar to the U = symmetric case 
considered in Sec. |^. This is qualitatively the case, and a very rapid narrowing of the 
central peak can be seen in the inset of the top panel in figure |^ for A > Ac where 
Xc = ^JUuJo/2 = 0.05. 
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Figure 5. Same as figure |l|, but witli e/ = —0.05. 



For small A, we see the typical charge excitations of the symmetric Anderson model 
at cij = ±f//2. For larger values of A, we see in general four high-energy peaks. The 
positions and relative weights of these peaks can be understood in terms of the diagram 
shown in figure |], where we plot the positions of the dominant excitation peaks taken 
from the atomic limit result displayed in the middle panel of figure |^. These are obtained 
from the atomic limit Green's function by replacing the multi-phonon peaks for each 
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Figure 6. Imaginary part of tlie phonon propagator for tire parameters of figure ^ 
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term in equation (|I5|) by a single peak with n excited phonons: 

(if + nujQ = ej + 2v from zero-occupied site (22a) 

ef — nuQ = ef from single-occupied site (226) 

e/ + f/ + nujQ = ef + U from single-occupied site (22 c) 

(Tf + U — nujQ = tf + U — 2v from double-occupied site (22 c?) 

The total weight in each peak is the sum of the weights of the contributing excitations 
in the atomic limit. The size of the symbols in each curve indicates the corresponding 
weight factor. One observes the transfer of weight with increasing A from the excitations 
in which the site is initially singly-occupied to those of the zero- and doubly-occupied 
states as U changes sign. In the approach to the regime where these levels cross, all 
four peaks can be clearly distinguished in figure 0. 

The imaginary part of the phonon Green's function is shown in figure ^ We see 
that while U > charge fluctuations are suppressed and there is no soft phonon mode 
ujx. This mode only emerges when U becomes negative. Its appearance clearly correlates 
with the collapse of the central peak. 

For comparison, the perturbational results corresponding to the lowest-order 
diagrams in U and A are shown in the lower panel of figure |^. It only gives a good 
approximation to the NRG results for A = 0, and there is clear disagreement for finite 
A as in the U = case discussed in Sec. El. 
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Figure 7. Spectral density for U = 0.1, e/ = -0.05, T/ = 0.1 ^ A « 0.016 and 
ujQ — 0.05. The respective values of A are given in the graph. The upper panel shows 
NRG results, the two lower panels the atomic limit and the perturbation theory (PT) 
results. For the atomic limit calculation, the height of the peaks is a measure of the 
respective spectral weight of the excitation. 



6. Asymmetric U Model 



In Figs. |T0| and |TT|, we plotted the spectral density for U = 0.2 and U - 
{ef = -0.05,1/ = 0.1 ^ A = 0.016, Uo = 0.05). In the first case. 



= 0.05, respectively 
for A = 0, we are 
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Figure 9. Imaginary part of the phonon propagator for the parameters of figure ^ 



in the Kondo regime and a narrow Kondo resonance at the Fermi level can clearly be 
distinguished from the atomic level below. 

For intermediate values of A, the spectrum again involves four dominant peaks 
whose relative positions and weights can be understood on the basis of the diagram on 
the left in figure 0. These excitations correspond exactly to the ones in equation (21). 
The transfer of weight with increasing A is rather different: in the intermediate regime, 
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Figure 10. Same as figure 0, but with U — 0.2. 



0.04 < A < 0.08, weight is transferee! from those excitations which involve an initial 
single-occupied state to those in which the initial state is unoccupied. As this transfer 
of weight develops, the Kondo resonance disappears as can be seen clearly in the inset 
of the top panel of figure llO . 



In the results with U = 0.05, figure |Tl| and the diagram on the right in figure 0, 
only one peak is to be seen for A = 0. reflecting the fact that for this set of parameters, 
we are in the mixed- valence regime. The change in high-energy features with increase 
of A can again be understood in terms of the transfer of weights between the dominant 
excitations as indicated in figure |T2|. For large A, the spectral weight accumulates for 
excitations from the double-occupied local level. Consequently, the peak moves to lower 
energies in a similar way to that shown in figure ^ for ?7 = (see discussion in Sec. 
The crossover occurs for smaller values of A compared to the case shown in figure ITD 
due to the smaller value of U. 
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Figure 11. Same as figure ^ but with U = 0.05. 



In both these cases, the charge fluctuations are small, and there is no evidence of a 
soft phonon mode developing in the phonon spectrum. Also, perturbation theory does 
not give reasonable results for either one of these cases. 

7. Conclusions 

We begin with an overview of the detailed results presented in sections 3-6, and look 
first of all at the results for the symmetric model. When A = and we move from 
U = to U = 0.1, we move from an intermediate valence regime, in which the spin 
and charge fluctuations are of the same order to the Kondo regime, in which the charge 
fluctuations are suppressed. The central resonance narrows considerably from a width 
2 A to one determined by the Kondo temperature Tk, and atomic-like high energy peaks 
begin to develop. When we switch on, and increase, the phonon coupling the central 
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Figure 12. Positions of the dominant peaks for the atomic limit calculation from 
Figs. |l^ and |ll|. The size of the circles shows the total weight for the respective type 
of excitations (see text). 
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Figure 13. Imaginary part of the dynamic charge (xc) and spin susceptibilities (xs) 
for the symmetric parameters of figure]^ {U = 0.1, e/ = —0.05, V — 0.1). 



peak initially broadens and then for large plionon coupling strength narrows markedly. 
The initial broadening can be understood as mainly due to the reduction in the effective 
electron interaction from U U = U — 2v, causing an increase in Tk, which depends 
exponentially on U. However, the atomic-like peaks at ~ ±U/2, seem to be largely 
unaffected. As U falls below the value ttA, the relative strength of the charge fluctuations 
increases, and we have a regime which is similar to the intermediate valence regime as 
the spin and charge fluctuations are of the same order, though high energy features 
associated with the atomic limit persist. Finally we move into the regime where U 
reaches zero and then becomes negative. In this limit we see that the spin fluctuations 
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are largely suppressed as pairing occurs and the impurity state is either doubly-occupied 
or empty. Excitations from these states are associated with broadened high energy 
features above and below the Fermi-level. The central resonance narrows exponentially. 
One can interpret the narrow resonance in the strong-coupling regime as analogue of the 
Kondo resonance in a negative-?/ model. In this regime the spin fluctuations have been 
largely eliminated but there remain long-lived charge fluctuations between the double 
and zero occupation states at the impurity site. The double and zero occupied states 
can be interpreted as the up and down configurations of a local isospin. In the a — t 
model an exchange coupling of such an isospin to the conduction electrons is invoked to 
mimic one of the channels in the two channel Kondo model [T^ . 

The reversal of roles of the spin and charge fluctuations can be seen in the results 
for the corresponding dynamic susceptibilities shown in figure |l^ for U = 0.1. For A = 
there is a low energy peak in the spin susceptibility that gradually gets suppressed for 
increasing A and moves to higher frequencies. The reverse trend is seen in the charge 
susceptibility where a low energy peak develops and becomes very narrow in the negative 
U range. 

We can also interpret the narrow central resonance as the zero-phonon excitation 
which carries very little weight in the strong-coupling regime because of the high 
probability of exciting a phonon when changing the impurity occupation. The narrow 
peak in the charge susceptibility is reflected in the phonon spectrum as low energy mode 
at ujx- 

Whereas the first of these interpretations would give a resonance that exponentially 
narrows as e~'^'/^'^^ with U — U ^ —\^/ujq^ we found that it narrows more rapidly than 
this expression would indicate (see figure The expression for the narrowing involving 
the soft phonon mode, ~ exp(— w/cua) with ux ~ exp(— A^/ciJq) gives a better fit to data. 

In the results for the asymmetric model the low energy phonon mode does not 
develop, and there is no narrow resonance at the Fermi-level other than the original 
Kondo peak in the case U = 0.2 at zero or small phonon coupling. In the large A regime 
{U < 0) in the two cases studied there is either a broad peak below the Fermi-level 
corresponding to double occupation, or one above the Fermi-level corresponding to zero 
occupation. For asymmetric conditions the double and zero occupied states are not 
equally likely, and the local level position e/ for e/ 7^ acts like a magnetic field on 
the isospin suppressing the 'Kondo' resonance. The same behaviour is observed in the 
negative-f/ Anderson model without phonons. 

Our conclusions are very much in line with those of Hewson and Newns in their 



study of the spinless model [|1^, |T4|. A direct comparison with their work is not 
possible because no f/-term can be induced by the phonons in the spinless case. Their 
conclusion was that a significant narrowing of the local resonance could only be induced 
by strong phonon coupling if the impurity level lies at, or very close to the Fermi-level 
(intermediate valence regime), when the zero or single occupation states of the local 
level have the same energy. This closely resembles the situation in the model with spin, 
but in terms of doubly-occupied and the empty states. 
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We have surveyed the physically interesting parameter regimes of the model using 
the NRG, and related our results to atomic limit calculations and perturbation theory. 
The relevant parameter regime for magnetic impurities systems is expected to be large 
U, with a phonon coupling such t> -C f/, and a;o < A. We have taken a relatively large 
value of uq {luq ~ ttA) in order to see the features in the spectrum clearly and save 
computational time (see Appendix B), but results for smaller values of loq were found 
to be qualitatively similar. 

Dynamical response functions for impurity systems are difficult to observe 
experimentally, so it would be interesting to calculate the effects of the phonon 
coupling on the thermodynamic behaviour which could be more readily compared with 
experiment. This we propose to do for a separate publication. As a more realistic model 
for magnetic impurities we also intend to include a direct coupling of the phonon mode 
to the hybridization term, because this term should be rather sensitive to local lattice 
vibrations. 

Having taken the first step of calculating the dynamic properties of the impurity 
model with phonons, we now propose to extend these calculations to lattice models 
within the framework of the dynamical mean-field theory. These calculations 
should be particularly relevant to high-Tc superconductors where recent photoemission 
experiments give evidence of strong electron-phonon coupling The electron- lattice 
coupling is also known to be important for the manganites [O]. The calculation could 



be extended to a model for these systems by including a coupling to a local spin |22 
which would require only a generalization of the impurity part of the Hamiltonian. The 
calculation could also be extended to include a dispersive phonon mode. This would 
require the introduction of a second chain of states, but then only a lower proportion of 
states could be kept in each iteration. 
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Appendix A. Equations of motion 

We can use the standard equation of motion for the Fourier transform of the double-time 
Green's function, 

{{AB))^ = {[A,B],) + {{[A,HUB))^ (A.l) 
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To derive the expression for the self-energy given in equation , we choose A = and 
B = fl and = + as these are fermion operators. The equations of motion give 

- es)G.{uj) = 1 + U{{fjlU, fl))^ + A((^(6t + 6), /]))^ + ^ V,((c^^, /t))^ (A.2) 

k 

(^-ek)((ck.,/])). = l^kG'.M (A.3) 

Hence 



where 



We define 



{uj-ef- /\{lo))G,{lo) = 1 + UF,{lo) + \M,{uj) (A.4) 
^.M = {{fJlf.. fDU M^ico) = iiUb^ + b), /t)). (A.5) 

so that S^(t<;) = TF^{uj) + S^(cj). 

For the boson Green's function D{uj) we take A = B = l)^ + h and rj = —. The 
equation of motions of motion give 

uDico) = LUo{{{-b^ + b), (6t + b)))^ (A.7) 
uiii-b^ + 6), (fet + b)))^ = 2 + 2A((g, (fet + 6))^ + ,;oD(cu) (A.8) 

Hence, 

D{uj) = D^iio) + AD°M((d, (6t + 6)))^ (A.9) 
where O = J2a flfa " 1- Taking the equation of motion for the right-hand operator, 

- uj{{d, (6t + b))))^ = uo{{d, + b)))^ (A.IO) 

- LoiiO, i-b^ + b)))^ = 2A((6, 6))^ + uo{{d, (fet + b)))^ (A.ll) 
((O, (6^ + b)))^ = D\u) + XD\u){{d, 6))^ (A.12) 

Hence the result, 

D{uj) = D\uj) + A2Z}°H((6, d))^D\uj) (A.13) 
Appendix B. Matrix elements and NRG calculation 

To evaluate the Green's function Mo-(u;), and hence the self-energy S^(co'), we need 
to calculate the matrix elements of the operator f„{b'^ + b). We can classify the states 
according to total charge and total spin quantum numbers Q,S and 5*2, and an additional 
index w or r, and we define reduced matrix elements using the Wigner-Eckart theorem 
in the usual way. The bose operators only enter explicitly in the calculation at the 
impurity site. Let the basis states used for this calculation be denoted by |Q, S", 5*2, n, r) 
where n denotes the occupation number for the bose states, and r is an additional 
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index. After diagonalization the exact eigenstates at the impurity site can be written in 
the form \Q, S, S^, w), where w no longer denotes the number of bosons, which are not 
conserved by the interaction term, so w is a general index. We need the reduced matrix 
elements between these eigenstates, {Q, S, w\ \f^{b^ + b)\\Q', S', w'). If 

\Q\S\w')=Y.UQ,,s'{w\n'r')\Q\S\n\r') (B.l) 

n' 

then 

{Q,S,w\\fAb^ + mQ',S',w') = 

= E UQ,s{w,n,r)UQ,Mw',n'r'){Q,S,n,r\\Ub^ + b)\\Q',S',n',r') 

n,n',r,r' 

= E UQ,s{w,n,r){^UQ,,s'{w',{n~iy)+ ^^'^^ 

n,r,r' 

+V^Uq>,s'{w', {n + iy)){Q,S\\fJ\Q',S') 

as {Q, S,n,r\\f^\\Q' , S' ,n,r') is independent of n. Hence we can deduce the values of 
{Q, S, w\\f^{b^ + b)\\Q', S', w') after diagonalization at the impurity site. Once we have 
these initial values the equations at all later steps will take the same form as those for 
{Q,S,wm\Q',S',w'). 

The spectral density Ccr{uj) of the Green's function Mo-(co') is given by 

CM = 1 + b^)f^\n'){n'\fl\n)6{u; - - E^)){e-^^- + e"^^-') (B.3) 

n,n' 

where \n) = \Qn, Sn, Sz,n,Wn) and Z is the partition function. When expressed in 
terms of reduced matrix elements it takes exactly the same form as the spectral density 
for Ga except that the reduced matrix element {Q, S,w\\f^{b^ + b)\\Q', S',w') replaces 
{Q,S,w\\a\Q', S',w'). 



NRG parameters: 

The NRG calculations were performed using the discretization parameter 1.7 < A < 2., 
and keeping at least 600 states in each iteration (see Ref. P, 0). The number of 
eigenstates of an impurity site coupled to a local phonon mode (the atomic limit (p!5|)), 
is already infinite. So a cutoff of higher-energy states is already necessary at the initial 
step of the NRG procedure. The probability distribution function to find n excited 
phonons in the system has a clear maximum around n = n = v/uq, as can be read off 
equation (p!5|), and falls off rapidly for larger values of n. So in the spirit of the NRG it 
is sufficient to keep only states with n < ncutofr where we chose an initial ncutofr = 4n. 

Appendix C. Derivation of equations and (|^) 

Here we give details of the explicit removal of the phonon terms in the Hamiltonian to 
derive equations (||) and (H). We will factorise the time-ordering operator T in equation 
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(0) so that T = TfTb, where Tp operates purely on the Fermi operators and Tb purely 
on the bose ones. As the interaction term with the phonons, 

i/B^(r) = A6(r)(6te-«" + 6e-^o"), (C.l) 

is linear in b and b\ we can apply the well-known Baker-Hausdorff formula to deduce 

We now take the expectation value of e~Jo -^i^t^'^)'^'^ with respect to the boson states 
using (e"^) = e2^^^\ applicable when A is a linear combination of b and b^ [17|. We find 
in the exponential, 

^ /' rd{r)d{r'){^{r-r')+^{r'-r))drdr'~ f HUT)dr, (C.3) 
/ Jo Jo Jq 

where H^^{t) is the purely fermion interaction term, and 

7(r - t') = n(cuo)e""("-"') + (1 + n{uJo))e-^"^^-^'\ (C.4) 
We also have that 

[/^int(r),i/int(r')] = ^6(r)6(r')(e--°(— ') - e^^^^-^'^), (C.5) 

subject to the fermion time-ordering within the complete expression. 
Finally, we find 

r/3 A, 



where 



-ujo\t-t'\ ^ujo\t-t'\ 



D\t -t')=[ + — ^ , (C.7) 



which is the imaginary time representation of the phonon Green's function. 
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